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Abstract 

A revised version of the quaternion approach for numerical integration of the 
equations of motion for rigid polyatomic molecules is proposed. The modified ap- 
proach is based on a formulation of the quaternion dynamics with constraints. This 
allows to resolve the rigidity problem rigorously using constraint forces. It is shown 
that the procedure for preservation of molecular rigidity can be realized particu- 
larly simply within the Verlet algorithm in velocity form. We demonstrate that 
the presented method leads to an improved numerical stability with respect to the 
usual quaternion rescaling scheme and it is roughly as good as the cumbersome 
atomic-constraint technique. 
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1 Introduction 



Many models of statistical mechanics deal with systems composed of classical rigid 
molecules. The method of molecular dynamics (MD) is widely applied for studying such 
systems. All known MD techniques appropriate to simulate molecular liquids can be split 
into three main approaches. In the first two approaches, time evolution of the system 
is considered in view of translational and rotational motions. These approaches differ 
between themselves by parameters which are used to represent the rotational degrees of 
freedom. In the classical scheme [1, 2], an orientation of the molecule is defined in terms 
of three Eulerian angles. As is well known [3] , the equations of motion are singular in this 
case. To avoid the singularities, Barojas et al [4, 5] have used two different sets of the 
Eulerian angles, each of which is applied in dependence on the orientation of the molecule. 
However, this procedure involves additional complex transformations with transcendental 
functions. 

In the second approach the rotation motion of a molecule is described without involving 
Eulerian angles. Cheung [6] has shown how to remove the singularities using special 
properties of diatomic molecules. For these molecules, Singer [7] has derived rotational 
equations of motion in terms of radius-vectors passed from one atom to another within 
the same molecule. An extension of this scheme to triatomic molecules was considered 
also [8, 9]. An alternative scheme has been proposed by Evans et al [10-12], where using 
so-called quaternions [3, 13, 14] leads to a singularity free algorithm for rigid polyatomics. 

In the third approach, proposed originally by Ryckaert et al [15], the cartesian equa- 
tions of motion are applied with respect to individual atoms. The total force on a particle 
appears as the sum of the force deriving from the potential energy and the force arising 
due to holonomic constraints. These atomic constraints must be in part rigid bonds in 
part linear relations to provide the rigidness of arbitrary polyatomics [16]. 

Apart from removing singularities a benefit derived from the last two approaches lies 
in the avoidance of time-consuming trigonometric functions. However, on integrating 
the equations of motion numerically, one additional difficulty appears here, namely, the 
problem of exact conservation of molecular rigidity. In the usual integration algorithms 
the rigidity can not be conserved with the precision better than that of evaluating the 
atom trajectories. For overcoming this drawback, it is necessary either to perform the 
rescaling of quaternions [11] or, within the atomic-constraint technique, to find solutions 
for a complete system of nonlinear equations [15, 16] at each step of the integration. 
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From the aforesaid, a natural question appears about the existence of a scheme which 
is free of all these drawbacks and yet has all advantages inherent in the mentioned above 
approaches. In the present paper we develop the idea of using quaternions to treat rota- 
tional motion. Section 2 is devoted to a general formulation of the quaternion dynamics 
with constraints. Applications of this approach to particular algorithms are considered 
in Sec. 3. The problem of how to adapt the Verlet algorithm to integrate the quaternion 
equations of motion is also solved there. In Sec. 4 various approaches are compared and 
discussed. Some concluding remarks are given in Sec. 5. 



2 Quaternion dynamics with constraints 

We consider a system of N identical rigid molecules with mass m, composed of M 
point atoms. In the molecular approach, evolution of the system in time is separated 
into translational and rotational motions. The translational motion is applied to the 
molecule as a whole and can be described by the 3N (i = 1, . . . , N) Newton equations 
= ^jw^ Fij(\ r i ~ r jl)> where and rf are the positions of the center of mass and 
atom a of molecule % in the laboratory fixed coordinate system L, respectively, and F'f- 
are the atom-atom forces between two different molecules. 

In order to analyze rotational motions, we introduce the sets e = (e 1 ,e 2 ,e 3 ) and 
(u\, u 2 , u l 3 ) = Aje of orthogonal unit vectors characterizing the L-system and the moving 
coordinate system S* attached to molecule i, respectively, where Aj is a rotational matrix. 
The angular velocity f2 { of the i-th molecule is defined as du l a / dt — £2iXu % a . The principal 
components, Q{u\ + Q l 2 u 2 + Q\u\ = f2 i: of angular velocities (i — 1, . . . , N) obey the 3N 
Euler equations [1]: 

Ja ^dT = K{t) + ( J/3 " J -y) Q W) Q W > (!) 

where (a, (3, 7) = (1,2,3); (2,3,1) and (3,1,2). Here J 1: J 2 and J 3 are the moments 
of inertia along principal axes of the molecule, Ylf-^ = ^i e i + ^2 e 2 + k l 3 e 3 = 

K\u\ + K\u % 2 + K\u\ is the torque exerted on molecule % with respect to its center of 
mass due to the interactions with the other molecules, K { = Ajfcj and 81 — r° — r^. Let 
A a = (A", Aj, A3) be a vector-column of positions for atom a within the molecule in the 
S l -system, i.e., 81 = Alu\ + A 2 u l 2 + A 3 u\. Then the positions of atoms in the L-system 
at time t are rf(t) = ri(t) + A^~(t)A a , where A + denotes the matrix transposed to A. 
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It is a common practice to define an orientation of the S*-system with respect to 
the laboratory frame in terms of three Eulerian angles. A numerical integration of the 
corresponding equations of motion has been performed in early investigations [1, 2]. As 
was soon realized, however, this procedure is very inefficient because of the singularities 
whenever the azimuthal angle of a molecule takes the value or n [11]. It has been shown 
in later investigations [10, 11] that at least four orientational parameters per molecule 
(quaternion) must be used to avoid the singularities. 

The orientational matrix Aj = A(q i ) in terms of the quaternion q i = Tfo, C»j X») i s 
given by [10, 11]: 



MQi) 



-2(&77i + CiXi) 
\ 2(^0 - ZiXi) 



%-vf-Q + x* 

+ ViXi) 



2(^0 + ZiXi) ^ 

-^f-vf + Cf + xl J 



(2) 



and the time derivative of the quaternion is expressed via principal components of angular 
velocity as follows 



Qi = 



Vi 
Ci 

V Xi ) 



Ci Xi 



Vi 6 

Xi Ci Ci Vi 

Ci Vi Xi Ci 

V -Vi & -0 Xi J 



( n{\ 

V o ) 



^Q(Qi)^i , 



(3) 



where the matrix Q is the function of q r It is worth to underline that the matrix A 
is a rotational one if the quaternion satisfies the equality qf = Q + rjf + Q + xj — 1- 
Differentiating the relation (3) over time yields 



Qi = \ (Q(4i)Oi + Q(<?JsV) , 



(4) 



where f2j = 2Q~ 1 {q, j )q i . It is trivial to find that the inverse matrix Q _1 = Q + since 
the 4x4 matrix Q is orthogonal when q\ — 1. We have augmented the angular velocity 
vector to involve square matrices using the result q^ q i = + ViVi + CiCi + XiXi — 
which follows from the equality q\ — 1. 

Then, using the Newton and Euler equations (1), we obtain the coupled set of 7N 
second-order differential equations of motion J^({rj, r i: q i: q^ q\}) = in terms of the 7N 
generalized coordinates {?%, q { }. If an initial state {vj(to)) ^i(^o)> Qi{to), Qi{to)} is specified, 
the evolution {ri(t), <?j(t)} of the system can be unambiguously determined. 
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Let us look for an analytical solution of the quaternion equations of motion by writing 

where qf\to) denotes the p-fold time derivative of q, L at time t . It is easy to check 
from the structure of equation (3) that arbitrary-order time derivatives of the quaternion 
constraint a^t) = qf(t) — 1 = are equal to zero, i.e., q^q^ = 0, q\ + q i -q i = and so on. 
Therefore, if all terms (P — > oo) of the Taylor's expansion (5) are taken into account and 
initially all the constraints are satisfied, <7j(£ ) = 0, they will be fulfilled at later times 
as well. In practice, however, the equations of motion are not solved exactly, so that 
these constraints will only be satisfied approximately. Let the integration algorithm used 
involves an error in the coordinates of order At p+1 , where At — t — to is the time step. 
In the simplest case of Taylor's expansion (5), this is the order of the first omitted term. 
Then the same order of uncertainties will be accumulated at each time step in conservation 
of the molecular rigidity, i.e., a^t) = 0(At p+1 ). In such a case, molecules are collapsed 
or even destroyed in time. In the usual version [11] of the quaternion method, to achieve 
the required rigidity at all times, it was proposed to multiply each quaternion component, 
associated with the same molecule, on the common factor 1 j \J~q? at every time step of 
the numerical integration (the so-called rescaling scheme). 

We consider now the question how to replace the crude quaternion renormalization 
by a more natural procedure in the framework of a systematic approach. The fact that 
quaternion components are not independent, requires, in general, the necessity of intro- 
ducing additional forces, namely, f^t) = —\ i (t)'V q .<Ji(t) = — 2X i (t)q i (t), which appear as 
a result of the constraints. These virtual quaternion-constraint forces should be added to 
the equations of motion (4) and, as a consequence, they modify the solution (5) as follows 

«.<*) = £ «i"<*o)*4f£ + £ f i r y o) (±-hl , (6) 

p=0 P- p=2 P- 

where /?" 2) (t ) = -^ElZlC^ 2 \\ k) (t )q^ 2 ' k) (t ) denote the (p - 2)-fold time deriva- 
tives of constraint forces, \\°\to) is a value of the Lagrange multiplier Aj(t), and \\ k \t ) 
are its fc-fold time derivatives (k — 1, . . . , P — 2) at time t . Differentiating (6), we obtain 
/-fold time derivatives (I — 0, . . . , P — 2) of q i at time t: 

?<^)=£^U)^^-2 £ £^A<%)<r 2 -^)^-^ . (7) 

p=l ^P >' p=max{2,Z} fc=0 vA* /" 
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In order to computer P — 1 unknowns (to), we have merely to exploit the information 
contained in the constraint (Ji(t) = Qi(t) 2 — 1=0. As this holds at any time, at least the 
first P — 2 time derivatives of cXj(t) must vanish. Then one obtains {p = 1, P — 2): 

= -^w =* P f:ci- 1 q?\t) q <r k \t) = o . (8) 

QI fc=0 

In view of explicit expressions (7), the conditions (8) together with the basic con- 
straints qf (t) = 1 constitute a system of P — 1 nonlinear equations per molecule with 
respect to the same number of unknowns X\ k \t ). The equations can be linearized and 
solved in a quite efficient way by iteration. This is justified for At — > because then 
the terms nonlinear in X[ k \t ) are small. Thus the iteration procedure can be initiated 
by substituting X[ k \t ) = in all nonlinear terms, and iterations always converge rapidly 
to the physical solutions X\ k \t ) ~ At p ~ k ~ l . The contributions of quaternion-constraint 
forces into the quaternion dynamics (6) are of order At p+1 , i.e., the same order as un- 
certainties of the integration algorithm (5), but the rigidity is now fulfilled perfectly for 
arbitrary times in future. It is worth emphasizing that these forces are imaginary and 
depend on details of the numerical integration in a characteristic way, contrary to the real 
bond forces in the atomic-constraint dynamics [15, 16]. They vanish if the equations of 
rotational motion are solved exactly. 



3 Applying actual algorithms 

3.1 Integration within the Gear method 

Usually, the Gear predictor-corrector algorithm [17, 18] is applied to integrate the 
equations of rotational motion [1, 11, 12]. In particular, it has been used [11, 12] for the in- 
tegration of the quaternion equations. Within the Gear method the quaternions and their 
time derivatives are predicted using the Pascal triangle qf\t + At) = J2p=i qf^ (t) fpliy. ? 
where / = 0, . . . , P and P is the order of the algorithm. Further, they are corrected one 
or more times, using new values of torques as well as rotational velocities and their time 
derivatives which are predicted and corrected simultaneously with quaternion variables. 

The Gear method can be modified within the quaternion-constraint dynamics as fol- 
lows. To simplify notations, we choose the fourth order scheme (P = 4) (the extension to 
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arbitrary orders is trivial). Let q\ (t + At) (as well as q\ (t)) be already defined quan- 
tities after the last step of the corrector procedure. Then, according to the constraint 
formalism, the variables q^t + At), q { (t + At) and q\(t + At) (/ = 0, 1, 2) transform into 

q\{t + At) = q t (t + At) + fAt)At 2 /2 + f^At'/Q + f^A^/U , 

q\{t + At) = q t (t + At) + f t (t)At + f^At 2 ^ + f t (t)At 3 /Q , (9) 

q\{t + At) = q\(t + At) + Ut) + f^At + f^At 2 ^ , 

where = -2A i q fi (t), U{t) = -2(^^(0 + ^(t)), = -2(A i q i (t) + 2A i g i (t) + 

Xiq^i)) and Aj, Aj, Aj are values of the Lagrange multiplier and its first and second time 
derivatives at time t. The expressions (9) present, in fact, (in somewhat other notations) 
a particular case (P = 4) of generalized equations (7). Therefore, the three unknowns Aj, 
Aj and \ are found solving by iteration the system of three nonlinear equations 

q'f = l, q'i-q'^0, q' 2 + q\ • q\ = . (10) 

As in the general case (8), the iteration procedure is initiated by putting Aj = Aj = Aj = 
in nonlinear terms, and unknown quantities quickly tend to the physical solutions Aj ~ 
At 3 , Aj ~ At 2 and Aj ~ At. 

3.2 Verlet algorithm in velocity form 

There are the well-known group of integrators comprising Verlet [19], leapfrog [20], 
velocity Verlet [21] and Beeman [22] methods. Due to their simplicity and exceptional 
numerical stability they play an important role in the classical methodology of molecu- 
lar dynamics. All or some of these methods are always described and compared in any 
modern textbook [13, 14, 20, 23, 24]. However, the mentioned above approaches, being 
constructed initially for the integration of Newton's equations for translational motion, 
are not necessarily applicable directly to rotational dynamics. To our knowledge, only the 
leapfrog method has its versions for rotational motion [13]. The reason of such a situation 
is that contrary to translational dynamics, the second time derivatives of variables, asso- 
ciated with rotational degrees of freedom, may depend on their first time derivatives. In 
our case the pattern is complicated additionally by the necessity of including constraints 
in the equations of motion. We shall show now how to solve these problems within the 
Verlet algorithm in velocity form. 
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Let {ri(t),ri(t), q^t), <7i(^)} be a spatially- velocity configuration of the system at time 
t and {(Ji(t) = q^t) 2 — 1 = 0, &i(t) = 2q t - q t = 0}. The translational part {r i {t),r i {t)') of 
variables is considered within the Verlet algorithm in the usual way [21, 24], whereas the 
rotational variables {q^t), can be evaluated as follows. Using the principal torques 

Ki(t), we define angular accelerations i?j(t) and, therefore, second time derivatives q\(t) 
(4) on the basis of equations (1) for rotational motion. Then, taking into account the 
constraint forces f^t) = — 2X i (t)q i (t) yields 

Qi (t + At) = qi (t) + Ut)At + q\(t)At 2 /2 + + 0(At 3 ) . (11) 

The Lagrange parameters Aj are defined from the constraint relations <7j(t + At) = 
q 2 (t + At) — 1 = which constitute a single quadratic equation per molecule with the 
following solutions 



\ _ 1 



1 - q 2 At 2 /2 T Jl- q\AV - q^At* - (q 2 - qf) At A /4 



(12) 



where the time derivatives of quaternions are taken at time t. As can be verified easily, 
only the first solution is in self-consistency with the integration scheme. In the limit of 
small time steps, this solution behaves as \n — > q i -q i At/2, i.e., f^t) ~ At. Therefore, 
the constraint forces contribute into the quaternion dynamics (11) terms of order At 3 , i.e., 
the same order as numerical errors of the used algorithm, but the rigidity of molecules is 
now fulfilled exactly, i.e., Gi{t + At) = 0. 

And now we consider how to perform the second step 

s(t + At) = s(t) + (s(t) + s(t + At)) At/2 + 0(At 3 ) (13) 

of the velocity Verlet method, where s denotes a spatial coordinate. There are no problems 
to pass this step in the case of translational motion, when s = Tj and s = V{ is the 
translational velocity. However, the difficulties immediately arise for rotational motion, 
because then the second time derivative s can depend explicitly not only on the spatial 
coordinate s, but on the generalized velocity s as well. For example, choosing s = q { , 
we obtain on the basis of equations of motion (1) and (4) that q\(t) = q^q^t) , q { (t)) . 
In view of (13) this leads to a very complicated system of four nonlinear equations per 
molecule with respect to four unknown components of the quaternion velocity q i (t + At). 
It is necessary to note that analogous problems appear at attempts to apply the leapfrog, 
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usual Verlet and Beeman methods for rotational motion (even much more difficult in the 
last two cases). 

An alternative has been found in a rotational motion version [13] of the leapfrog 
algorithm. It has been suggested to associate the quantity s with the angular momentum 
Zj = A.f Li of the molecule in the laboratory system of coordinates, i.e., s = l i: where 
Li = ( J\Q\, -h^D = JS^i and J is the diagonal matrix of principal moments of 
inertia. Then the equation (13) is simplified, 

li(t + At) = k(t) + (ki(t) + ki(t + At)) At/2 + 0(At 3 ) (14) 

and, therefore, li(t + At) are easily evaluated using the torques ki(t + At) in the new 
spatial configuration {q^t + At)}. At the same time, new values for principal angular and 
quaternion velocities are obtained (when they are needed) using the relations f2i(t+At) = 
r l Ai(t + At)li(t + At) and q^t + At) = |Q(^(i + At))fli(t + At). 

Finally, we consider the third version of the velocity Verlet method for rotational 
motion. The idea consists in using angular velocities as independent parameters for 
describing the sate of the system in phase space. Then choosing s = fi^ and taking into 
account Euler equations (1), we obtain from (13) the following result 



(r a {n \t + At) = a i a (t)+ At 



K(t) + K(t + At) 



2J a 

+ - J 7 )(^(t)^(t) + Qi {n ' l \t + At)wl n ~ l \t + At 



(15) 



Unless Ji = J 2 = J3, the equations (15) are, in fact, the system of three quadratic 
equations per molecule with respect to the three unknowns f2 l a (t + At). The system (15) 
is relatively simple and can be solved by iteration (n = 1, 2, . . .) with O^ ' \t+At) = Q^it) 
as an approximation of zero order. A few iterations is sufficient for actual time steps to 
find the desired solutions with a great precision. 

From a mathematical point of view, all the three representations s = q^U or 17 j 
are completely equivalent, because the knowledge of an arbitrary quantity from the set 
(<ji, h, fi/) allows us to determine uniquely the rest of two ones. In the case of numerical 
integration the pattern is qualitatively different, because the investigated quantities are 
evaluated approximately. The choice s = q { can not be recommended for calculations due 
to its complexity. Moreover, it has yet a major disadvantage that the equality q { - q t = 
is broken at time t + At, whereas this equality remains in law by construction in other 
two cases. 
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The case s = li is the most attractive in view of the avoidance of nonlinear equations. 
Computations show (see Sec. 4), however, that the best numerical stability with respect 
to the total energy conservation exhibits the third version, when s = The reason of 
this can be found, taking into account that a kinetic part, \ Ya=\{J\{^\) 2 + ■h{Q % 2) 1 + 
J 3 (i?3) 2 ), of the total energy is calculated directly from principal angular velocities. At 
the same time, to evaluate angular velocities within the previous approach the additional 
transformations fii = J -1 AjZj with approximately computed matrices Aj = A(qr i ) and 
angular momenta li are necessary. They contribute an additional portion into accumulated 
errors at the calculations of the total energy. 

4 Numerical results and discussion 

We now test our integration approach on the basis of MD simulations for a TIP4P 
model [25] of water. The simulations were performed in the microcanonical ensemble at 
a density of 1 g/cm 3 and at a temperature of 298 K. We considered a system of N = 256 
molecules in the reaction field geometry [26]. All runs were started from an identical 
well equilibrated configuration. Numerical stability was identified in terms of the relative 
fluctuations S t = \J {{E — (E) t ) 2 ) t j (E) t of the total energy of the system during time t. 

Samples of the function $ t for the fourth-order Gear algorithm are presented in fig. 1. 
It can be seen easily from the figure that within the modified version, where quaternion- 
constraint transformations (9) are applied, the conservation of the total energy is im- 
proved considerably in comparison with that obtained within the original version of the 
algorithm, when no additional corrections are used. The system of nonlinear equations 
(10) were solved with the relative iteration precision of 10 -12 and the number of iterations 
never exceeded 5. At the same time, applying the rescaling scheme leads even to worse 
results than those given by the original version. Despite the fact that the Gear algorithm 
integrates the equations of motion very well at At < 1 fs, it has a very small region of 
stability with respect to time steps [16] and can not be used for At > 1.5 fs (see fig. Id). 
Moreover, the scheme predictor-corrector-corrector, which was chosen by us at At < 1 fs 
to provide an optimal performance, takes computational time per step twice larger than 
the Verlet integrator (in the case At = 1.5 fs, three corrector steps were used). 

As the atomic-constraint technique [15, 16] is intensively exploited and its advantages 
with respect to numerical stability are generally recognized, we have made comparative 
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test carrying out explicit MD runs using this method as well as the angular- velocity Verlet 
algorithm (15) within the quaternion-constraint dynamics (11). The usual value of the 
time step in simulating such a system is At = 2 fs [27]. These two approaches required 
almost the same computer time per step (96% being spent to evaluate pair interactions). 
The corresponding functions S t are shown in fig. 2 at four fixed values At — 1, 2, 3 and 4 fs. 
For the purpose of comparison the results obtained within the angular-momentum version 
(14) are also included in this figure (insets (a), (b)). The rescaling scheme (instead of the 
quaternion-constraint dynamics) within the angular- velocity Verlet version was considered 
as well. 

At the smallest time step (fig. 2a), all the approaches exhibit similar equivalence in 
the energy conservation. The time step At — 1 fs is too small, however, and impracti- 
cal in simulations because it requires too much computation time to cover the sufficient 
phase space. For larger time steps (fig. 2b-d), the total energy fluctuations increase drasti- 
cally with increasing the length of the runs within the usual quaternion rescaling scheme, 
whereas the atomic- and quaternion-constraint methods conserve the energy approxi- 
mately with the same accuracy. As far as the rescaling of quaternions has been chosen, 
it must be applied after each time step to achieve an optimal performance within the 
Verlet integrator. For example, if no corrections are performed, the rigidity of molecules 
is broken catastrophically (see the dotted curve in fig. 2b). As we can see, the modified 
quaternion approach always leads to improved results. Finally, we show in fig. 2a,b that 
the angular-momentum version of the Verlet algorithm is much more unstable and can 
be used at small (At < 1 fs) time steps only. For these reasons, taking into account also 
comments on the Gear algorithm, the crude renormalization procedure is generally not 
recommended and preference must be given to the modified quaternion approach within 
the angular-velocity Verlet integrator. Quite a few iterations (the mean number of itera- 
tions per molecule varied from 2 to 4 at At — 1 - 4 fs) was sufficient to obtain solutions 
to the system of nonlinear equations (15) with the relative iteration precision of 1CT 12 . 
This required a negligible small computation time additionally to the total time. 

The calculations have shown that the same level £ = 0.025% of energy conservation 
can be provided by the time steps of 2.1, 3.7 and 4.0 fs within the quaternion rescaling, 
atomic- and quaternion-constraint schemes. Therefore, the last two approaches allow a 
time step approximately twice larger than the quaternion rescaling method. A reason 
of this gain in time can be explained by the fact that the rescaling of quaternions is an 
artificial procedure. It involves an unpredictable discrepancy in the calculation of trajec- 
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tories of atoms at each step of the integration process and leads to significant deviations 
of the total energy. At the same time, the atomic- and quaternion(molecular)-constraint 
techniques provide the rigidity of molecules in a more natural way. 

5 Conclusion 

An alternative scheme to overcome the difficulties in simulations of rigid polyatomics 
has been proposed. The scheme uses the constraint formalism for treating quaternion dy- 
namics. As a result, the rigidity problem has been rigorously resolved, using quaternion- 
constraint forces. Although this introduces some extra transformations, but presents no 
numerical difficulties. In a particular case of the velocity Verlet algorithm, the constraint 
version of the quaternion approach allows one to perfectly fulfil the rigidity of molecules at 
each step of the trajectory without any additional efforts and loss of precision. Avoidance 
of the necessity to solve complex nonlinear equations for maintaining molecular rigid- 
ity should be a benefit of the presented approach with respect to the atomic-constraint 
dynamics. 

It has been corroborated by actual MD simulations that the quaternion rescaling 
method is much less efficient than the atomic- and quaternion-constraint techniques. The 
last both schemes seem to be comparable in efficiency The advantage of the modified 
quaternion approach is that it improves the energy conservation considerably at a little 
computational cost. 
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Figure captions 



Fig. 1. The relative total energy fluctuations as functions of the length of the sim- 
ulations performed within three versions of the Gear algorithm at four fixed time steps. 
Note, that all three curves are indistinguishable in (a). 

Fig. 2. The relative total energy fluctuations as functions of the length of the sim- 
ulations performed within the atomic-constraint technique and various versions of the 
velocity Verlet algorithm at four fixed time steps. The dotted curve in (b) corresponds 
to the usual quaternion approach without any additional corrections. Note, that three of 
four curves are indistinguishable in (a). 
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